30-09-2021_AdvancedData_ModSelection
Connor Hughes
9/30/2021
Load in packages…
Add in our data and make “Year” a catagorical predictor
#Read in Raw Data
TurtleData <- read.csv("Data/TidyNests.csv")
TurtleData$Year <- as.factor(TurtleData$Year)
TurtleData$X <- NULL
Make some Histograms for Nest Abundances and False Crawl Abundances
#Visualize our Counts Data. We can see that its not normally distributed, thats okay because we knew it wouldnt be
#Count Data falls under the Poisson distribution, so this visual makes sense to see
hist(TurtleData$NestAbundance)
### We'll Use This Information LaterTest our predictor variables for correlation among scales
Mean.Dep <- TurtleData[, c("MEAN.Depth.250", "MEAN.Depth.500",
"MEAN.Depth.750", "MEAN.Depth.1000")]
ggpairs(Mean.Dep)
SD.Dep <- TurtleData[, c("STD.Depth.250", "STD.Depth.500",
"STD.Depth.750", "STD.Depth.1000")]
ggpairs(SD.Dep)
Mean.Slp <- TurtleData[, c("MEAN.Slope.ON", "MEAN.Slope.250", "MEAN.Slope.500",
"MEAN.Slope.750", "MEAN.Slope.1000")]
ggpairs(Mean.Slp)
SD.Slp <- TurtleData[, c("STD.Slope.ON", "STD.Slope.250", "STD.Slope.500",
"STD.Slope.750", "STD.Slope.1000")]
ggpairs(SD.Slp)
Dists <- TurtleData[, c("Distance.to.0m", "Distance.to.10m",
"Distance.to.20m", "Distance.to.30m")]
ggpairs(Dists)
Unsuprisingly, we see corrilation among scales here for most of our data. This makes sense thinking about what we know about the ocean, and we do not want to use the same variable more than once anyway. So before we actually build our global model, we want to use AICc to see what scales these predictor variables are acting on our response variables. Lets do some scale selection! First lets build a bunch of mono-variable models.
options(na.action = "na.fail")
for (i in (2:25)) {TurtleData[,i] <- c(scale(TurtleData[,i], center = TRUE, scale = TRUE))}
NullModel <- glmer(NestAbundance ~ (1|Year), data = TurtleData, family = poisson)
DepthMean250 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.250, data = TurtleData, family = poisson)
DepthSTD250 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.250, data = TurtleData, family = poisson)
SlopeMean250 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.250, data = TurtleData, family = poisson)
SlopeSTD250 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.250, data = TurtleData, family = poisson)
DepthMean500 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.500, data = TurtleData, family = poisson)
DepthSTD500 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.500, data = TurtleData, family = poisson)
SlopeMean500 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.500, data = TurtleData, family = poisson)
SlopeSTD500 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.500, data = TurtleData, family = poisson)
DepthMean750 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.750, data = TurtleData, family = poisson)
DepthSTD750 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.750, data = TurtleData, family = poisson)
SlopeMean750 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.750, data = TurtleData, family = poisson)
SlopeSTD750 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.750, data = TurtleData, family = poisson)
Rugosity250 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.250 * STD.Slope.250, data = TurtleData, family = poisson)
Rugosity500 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.500 * STD.Slope.500, data = TurtleData, family = poisson)
Rugosity750 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.750 * STD.Slope.750, data = TurtleData, family = poisson)
Rugosity1000 <-glmer(NestAbundance ~ (1|Year) + STD.Depth.1000 * STD.Slope.1000, data = TurtleData, family = poisson)
DepthMean1000 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.1000, data = TurtleData, family = poisson)
DepthSTD1000 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.1000, data = TurtleData, family = poisson)
SlopeMean1000 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.1000, data = TurtleData, family = poisson)
SlopeSTD1000 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.1000, data = TurtleData, family = poisson)
Dist.to.0 <- glmer(NestAbundance ~ (1|Year) + Distance.to.0m, data = TurtleData, family = poisson)
Dist.to.10 <- glmer(NestAbundance ~ (1|Year) + Distance.to.10m, data = TurtleData, family = poisson)
Dist.to.20 <- glmer(NestAbundance ~ (1|Year) + Distance.to.20m, data = TurtleData, family = poisson)
Dist.to.30 <- glmer(NestAbundance ~ (1|Year) + Distance.to.30m, data = TurtleData, family = poisson)Then, lets pit all the scale models at each variable against one another.
depth.mean.aic <- list("DepthMean250" = DepthMean250, "DepthMean500" = DepthMean500,
"DepthMean750" = DepthMean750, "DepthMean1000" = DepthMean1000, "Null" = NullModel)
model.sel(depth.mean.aic)## Model selection table
## (Int) MEA.Dpt.250 MEA.Dpt.500 MEA.Dpt.750 MEA.Dpt.1000
## DepthMean1000 1.463 -0.2084
## DepthMean750 1.472 -0.1631
## DepthMean500 1.479 -0.107
## Null 1.485
## DepthMean250 1.485 -0.006744
## family df logLik AICc delta weight
## DepthMean1000 poisson(log) 3 -2914.479 5835.0 0.00 1
## DepthMean750 poisson(log) 3 -2955.517 5917.1 82.08 0
## DepthMean500 poisson(log) 3 -2993.886 5993.8 158.81 0
## Null poisson(log) 2 -3023.143 6050.3 215.32 0
## DepthMean250 poisson(log) 3 -3023.028 6052.1 217.10 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(DepthMean1000)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: NestAbundance ~ (1 | Year) + MEAN.Depth.1000
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 5835.0 5849.9 -2914.5 5829.0 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7972 -1.1773 -0.2841 0.7750 7.0389
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.1186 0.3444
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.46336 0.10995 13.31 <2e-16 ***
## MEAN.Depth.1000 -0.20845 0.01421 -14.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## MEAN.D.1000 0.027car::Anova(DepthMean1000)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: NestAbundance
## Chisq Df Pr(>Chisq)
## MEAN.Depth.1000 215.1 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
slope.mean.aic <- list("SlopeMean250" = SlopeMean250, "SlopeMean500" = SlopeMean500,
"SlopeMean750" = SlopeMean750, "SlopeMean1000" = SlopeMean1000, "Null" = NullModel)
model.sel(slope.mean.aic)## Model selection table
## (Int) MEA.Slp.250 MEA.Slp.500 MEA.Slp.750 MEA.Slp.1000
## SlopeMean250 1.457 -0.2367
## SlopeMean500 1.465 -0.2045
## SlopeMean750 1.463 -0.2291
## SlopeMean1000 1.466 -0.221
## Null 1.485
## family df logLik AICc delta weight
## SlopeMean250 poisson(log) 3 -2885.812 5777.6 0.00 1
## SlopeMean500 poisson(log) 3 -2927.465 5861.0 83.31 0
## SlopeMean750 poisson(log) 3 -2929.919 5865.9 88.22 0
## SlopeMean1000 poisson(log) 3 -2943.952 5893.9 116.28 0
## Null poisson(log) 2 -3023.143 6050.3 272.65 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(SlopeMean250)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: NestAbundance ~ (1 | Year) + MEAN.Slope.250
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 5777.6 5792.6 -2885.8 5771.6 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1863 -1.1104 -0.2273 0.8743 5.6166
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.1186 0.3444
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.45739 0.10997 13.25 <2e-16 ***
## MEAN.Slope.250 -0.23669 0.01442 -16.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## MEAN.Sl.250 0.030car::Anova(SlopeMean250)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: NestAbundance
## Chisq Df Pr(>Chisq)
## MEAN.Slope.250 269.5 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
rugosity.aic <- list("Rugosity250" = Rugosity250, "Rugosity500" = Rugosity500,
"Rugosity750" = Rugosity750, "Rugosity1000" = Rugosity1000, "Null" = NullModel)
model.sel(rugosity.aic)## Model selection table
## (Int) STD.Dpt.250 STD.Slp.250 STD.Dpt.250:STD.Slp.250 STD.Dpt.500
## Rugosity500 1.542 0.01304
## Rugosity1000 1.502
## Rugosity750 1.485
## Rugosity250 1.398 0.0298 -0.2719 0.07936
## Null 1.485
## STD.Slp.500 STD.Dpt.500:STD.Slp.500 STD.Dpt.750 STD.Slp.750
## Rugosity500 -0.2425 -0.1725
## Rugosity1000
## Rugosity750 -0.2771 0.265
## Rugosity250
## Null
## STD.Dpt.750:STD.Slp.750 STD.Dpt.1000 STD.Slp.1000
## Rugosity500
## Rugosity1000 -0.3159 0.3908
## Rugosity750 -0.06816
## Rugosity250
## Null
## STD.Dpt.1000:STD.Slp.1000 family df logLik AICc delta
## Rugosity500 poisson(log) 5 -2801.814 5613.7 0.00
## Rugosity1000 -0.0898 poisson(log) 5 -2886.044 5782.1 168.46
## Rugosity750 poisson(log) 5 -2888.796 5787.6 173.96
## Rugosity250 poisson(log) 5 -2917.066 5844.2 230.50
## Null poisson(log) 2 -3023.143 6050.3 436.61
## weight
## Rugosity500 1
## Rugosity1000 0
## Rugosity750 0
## Rugosity250 0
## Null 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(Rugosity500)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: NestAbundance ~ (1 | Year) + STD.Depth.500 * STD.Slope.500
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 5613.6 5638.6 -2801.8 5603.6 1075
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0761 -1.0635 -0.2950 0.7714 6.0333
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.1186 0.3444
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.54203 0.11047 13.959 < 2e-16 ***
## STD.Depth.500 0.01304 0.02375 0.549 0.583
## STD.Slope.500 -0.24246 0.03239 -7.486 7.08e-14 ***
## STD.Depth.500:STD.Slope.500 -0.17249 0.01903 -9.064 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) STD.Dp.500 STD.S.
## STD.Dpt.500 -0.014
## STD.Slp.500 0.068 -0.663
## STD.D.500:S -0.088 0.116 -0.384car::Anova(Rugosity500)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: NestAbundance
## Chisq Df Pr(>Chisq)
## STD.Depth.500 2.604 1 0.1066
## STD.Slope.500 141.213 1 <2e-16 ***
## STD.Depth.500:STD.Slope.500 82.154 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
dist.to.dep.aic <- list("Dist.to.0" = Dist.to.0,
"Dist.to.10" = Dist.to.10,
"Dist.to.20" = Dist.to.20,
"Dist.to.30" = Dist.to.30)
model.sel(dist.to.dep.aic)## Model selection table
## (Int) Dst.to.0m Dst.to.10m Dst.to.20m Dst.to.30m family df
## Dist.to.20 1.429 0.3331 poisson(log) 3
## Dist.to.30 1.433 0.3297 poisson(log) 3
## Dist.to.0 1.439 -0.3058 poisson(log) 3
## Dist.to.10 1.485 0.01826 poisson(log) 3
## logLik AICc delta weight
## Dist.to.20 -2737.011 5480.0 0.00 1
## Dist.to.30 -2774.551 5555.1 75.08 0
## Dist.to.0 -2793.227 5592.5 112.43 0
## Dist.to.10 -3022.285 6050.6 570.55 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(Dist.to.20)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: NestAbundance ~ (1 | Year) + Distance.to.20m
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 5480 5495 -2737 5474 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9521 -1.0876 -0.2367 0.7884 5.3526
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.1186 0.3444
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.42862 0.11002 12.98 <2e-16 ***
## Distance.to.20m 0.33306 0.01398 23.83 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Dstnc.t.20m -0.043car::Anova(Dist.to.20)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: NestAbundance
## Chisq Df Pr(>Chisq)
## Distance.to.20m 567.73 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now that we have selected the most relevant scale for each variable, build a new data frame and global model out of those variables at the most relevant scales to our response.
Turt.nest.data1 <- data.frame(nests = TurtleData$NestAbundance,
Year = as.numeric(TurtleData$Year),
muDep = abs(TurtleData$MEAN.Depth.1000),
muSlop = TurtleData$MEAN.Slope.250,
Rugosity = TurtleData$STD.Slope.500,
dist20 = TurtleData$Distance.to.20m,
shelfW = (TurtleData$Distance.to.20m - TurtleData$Distance.to.10m)
)
ggpairs(Turt.nest.data1, columns = 3:7) + theme_bw()
GlobalModel <- glmer(nests ~ muDep + muSlop + Rugosity +
dist20 + shelfW + (1|Year), data = Turt.nest.data1, family = poisson)
performance::check_collinearity(GlobalModel)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## muDep 1.15 1.07 0.87
## Rugosity 1.33 1.16 0.75
## shelfW 3.23 1.80 0.31
##
## Moderate Correlation
##
## Term VIF Increased SE Tolerance
## muSlop 5.14 2.27 0.19
## dist20 7.40 2.72 0.14GlobalModel_Dist <- glmer(nests ~ muDep + muSlop + Rugosity +
dist20 + (1|Year), data = Turt.nest.data1, family = poisson)
performance::check_collinearity(GlobalModel_Dist)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## muDep 1.03 1.02 0.97
## muSlop 2.41 1.55 0.41
## Rugosity 1.30 1.14 0.77
## dist20 2.79 1.67 0.36GlobalModel_Shelf <- glmer(nests ~ muDep + muSlop + Rugosity +
shelfW + (1|Year), data = Turt.nest.data1, family = poisson)
performance::check_collinearity(GlobalModel_Shelf)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## muDep 1.08 1.04 0.93
## muSlop 1.08 1.04 0.92
## Rugosity 1.32 1.15 0.76
## shelfW 1.22 1.11 0.82For now for simplicities sake, we can dredge our global model (while ignoring our multicoliniar variables)
dredge(GlobalModel_Shelf)## Fixed term is "(Intercept)"## Global model call: glmer(formula = nests ~ muDep + muSlop + Rugosity + shelfW +
## (1 | Year), data = Turt.nest.data1, family = poisson)
## ---
## Model selection table
## (Intrc) muDep muSlp Rgsty shlfW df logLik AICc delta weight
## 16 1.550 -0.1726 -0.1948 -0.1862 0.1990 6 -2638.878 5289.8 0.00 1
## 15 1.406 -0.2050 -0.1560 0.2149 5 -2654.294 5318.6 28.81 0
## 12 1.525 -0.1278 -0.2270 0.2547 5 -2670.354 5350.8 60.93 0
## 11 1.416 -0.2312 0.2612 4 -2679.250 5366.5 76.70 0
## 14 1.614 -0.2325 -0.2644 0.1869 5 -2717.685 5445.4 155.59 0
## 8 1.658 -0.2773 -0.1642 -0.3159 5 -2731.240 5472.5 182.70 0
## 13 1.421 -0.2248 0.2087 4 -2745.887 5499.8 209.98 0
## 7 1.429 -0.1762 -0.2786 4 -2773.047 5554.1 264.30 0
## 10 1.589 -0.1780 0.2782 4 -2778.305 5564.6 274.81 0
## 6 1.703 -0.3204 -0.3868 4 -2791.655 5591.3 301.51 0
## 9 1.438 0.2866 3 -2795.565 5597.2 307.32 0
## 5 1.440 -0.3463 3 -2847.544 5701.1 411.28 0
## 4 1.662 -0.2414 -0.2320 4 -2852.347 5712.7 422.90 0
## 3 1.457 -0.2367 3 -2885.812 5777.6 487.81 0
## 2 1.729 -0.2880 3 -2976.166 5958.4 668.52 0
## 1 1.485 2 -3023.143 6050.3 760.46 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'shelf.res <- dredge(GlobalModel_Shelf)## Fixed term is "(Intercept)"importance(shelf.res)## muSlop shelfW Rugosity muDep
## Sum of weights: 1 1 1 1
## N containing models: 8 8 8 8car::Anova(GlobalModel_Shelf)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: nests
## Chisq Df Pr(>Chisq)
## muDep 30.729 1 2.967e-08 ***
## muSlop 153.033 1 < 2.2e-16 ***
## Rugosity 56.643 1 5.227e-14 ***
## shelfW 176.732 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1And finally, lets plot our pairwise variables for our best model
ggplot(Turt.nest.data1, aes(muDep, nests)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.nest.data1, aes(muSlop, nests)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.nest.data1, aes(Rugosity, nests)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.nest.data1, aes(shelfW, nests)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
What does this look like biologically? 
LETS DO IT ALL AGAIN, this time with False Crawls
#Read in Raw Data
TurtleData <- read.csv("Data/TidyFC.csv")
TurtleData$Year <- as.factor(TurtleData$Year)
TurtleData$X <- NULLfor (i in (2:25)) {TurtleData[,i] <- c(scale(TurtleData[,i], center = TRUE, scale = TRUE))}
NullModel <- glmer(FCAbundance ~ (1|Year), data = TurtleData, family = poisson)
DepthMean250 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.250, data = TurtleData, family = poisson)
DepthSTD250 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.250, data = TurtleData, family = poisson)
SlopeMean250 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.250, data = TurtleData, family = poisson)
SlopeSTD250 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.250, data = TurtleData, family = poisson)
DepthMean500 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.500, data = TurtleData, family = poisson)
DepthSTD500 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.500, data = TurtleData, family = poisson)
SlopeMean500 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.500, data = TurtleData, family = poisson)
SlopeSTD500 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.500, data = TurtleData, family = poisson)
DepthMean750 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.750, data = TurtleData, family = poisson)
DepthSTD750 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.750, data = TurtleData, family = poisson)
SlopeMean750 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.750, data = TurtleData, family = poisson)
SlopeSTD750 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.750, data = TurtleData, family = poisson)
Rugosity250 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.250 * STD.Slope.250, data = TurtleData, family = poisson)
Rugosity500 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.500 * STD.Slope.500, data = TurtleData, family = poisson)
Rugosity750 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.750 * STD.Slope.750, data = TurtleData, family = poisson)
Rugosity1000 <-glmer(FCAbundance ~ (1|Year) + STD.Depth.1000 * STD.Slope.1000, data = TurtleData, family = poisson)
DepthMean1000 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.1000, data = TurtleData, family = poisson)
DepthSTD1000 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.1000, data = TurtleData, family = poisson)
SlopeMean1000 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.1000, data = TurtleData, family = poisson)
SlopeSTD1000 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.1000, data = TurtleData, family = poisson)
Dist.to.0 <- glmer(FCAbundance ~ (1|Year) + Distance.to.0m, data = TurtleData, family = poisson)
Dist.to.10 <- glmer(FCAbundance ~ (1|Year) + Distance.to.10m, data = TurtleData, family = poisson)
Dist.to.20 <- glmer(FCAbundance ~ (1|Year) + Distance.to.20m, data = TurtleData, family = poisson)
Dist.to.30 <- glmer(FCAbundance ~ (1|Year) + Distance.to.30m, data = TurtleData, family = poisson)
#Initial AIC Scale Selectionsummary(NullModel)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: FCAbundance ~ (1 | Year)
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 4975.6 4985.6 -2485.8 4971.6 1078
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4331 -1.0344 -0.2850 0.7805 6.3098
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.2325 0.4822
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1580 0.1536 7.541 4.67e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1depth.mean.aic <- list("DepthMean250" = DepthMean250, "DepthMean500" = DepthMean500,
"DepthMean750" = DepthMean750, "DepthMean1000" = DepthMean1000, "Null" = NullModel)
model.sel(depth.mean.aic)## Model selection table
## (Int) MEA.Dpt.250 MEA.Dpt.500 MEA.Dpt.750 MEA.Dpt.1000
## DepthMean1000 1.151 -0.1162
## DepthMean750 1.154 -0.09448
## DepthMean500 1.156 -0.06834
## Null 1.158
## DepthMean250 1.158 -0.01142
## family df logLik AICc delta weight
## DepthMean1000 poisson(log) 3 -2459.981 4926.0 0.00 1
## DepthMean750 poisson(log) 3 -2468.583 4943.2 17.20 0
## DepthMean500 poisson(log) 3 -2476.781 4959.6 33.60 0
## Null poisson(log) 2 -2485.816 4975.6 49.66 0
## DepthMean250 poisson(log) 3 -2485.567 4977.2 51.17 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(DepthMean1000)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: FCAbundance ~ (1 | Year) + MEAN.Depth.1000
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 4926.0 4940.9 -2460.0 4920.0 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4152 -1.0250 -0.2146 0.6770 6.5444
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.2325 0.4822
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.15124 0.15357 7.497 6.55e-14 ***
## MEAN.Depth.1000 -0.11616 0.01617 -7.182 6.86e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## MEAN.D.1000 0.012slope.mean.aic <- list("SlopeMean250" = SlopeMean250, "SlopeMean500" = SlopeMean500,
"SlopeMean750" = SlopeMean750, "SlopeMean1000" = SlopeMean1000, "Null" = NullModel)
model.sel(slope.mean.aic)## Model selection table
## (Int) MEA.Slp.250 MEA.Slp.500 MEA.Slp.750 MEA.Slp.1000
## SlopeMean250 1.141 -0.1845
## SlopeMean1000 1.148 -0.1527
## SlopeMean500 1.149 -0.1328
## SlopeMean750 1.149 -0.144
## Null 1.158
## family df logLik AICc delta weight
## SlopeMean250 poisson(log) 3 -2421.862 4849.7 0.00 1
## SlopeMean1000 poisson(log) 3 -2453.637 4913.3 63.55 0
## SlopeMean500 poisson(log) 3 -2454.050 4914.1 64.38 0
## SlopeMean750 poisson(log) 3 -2454.537 4915.1 65.35 0
## Null poisson(log) 2 -2485.816 4975.6 125.90 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(SlopeMean250)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: FCAbundance ~ (1 | Year) + MEAN.Slope.250
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 4849.7 4864.7 -2421.9 4843.7 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3620 -1.0238 -0.1832 0.7159 5.5979
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.2325 0.4822
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.14113 0.15358 7.43 1.09e-13 ***
## MEAN.Slope.250 -0.18451 0.01642 -11.24 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## MEAN.Sl.250 0.019rugosity.aic <- list("Rugosity250" = Rugosity250, "Rugosity500" = Rugosity500,
"Rugosity750" = Rugosity750, "Rugosity1000" = Rugosity1000, "Null" = NullModel)
model.sel(rugosity.aic)## Model selection table
## (Int) STD.Dpt.250 STD.Slp.250 STD.Dpt.250:STD.Slp.250 STD.Dpt.500
## Rugosity500 1.158 0.1239
## Rugosity1000 1.181
## Rugosity250 1.060 0.04737 -0.248 0.1006
## Rugosity750 1.168
## Null 1.158
## STD.Slp.500 STD.Dpt.500:STD.Slp.500 STD.Dpt.750 STD.Slp.750
## Rugosity500 -0.3455 -0.05202
## Rugosity1000
## Rugosity250
## Rugosity750 -0.1722 0.1544
## Null
## STD.Dpt.750:STD.Slp.750 STD.Dpt.1000 STD.Slp.1000
## Rugosity500
## Rugosity1000 -0.2266 0.2719
## Rugosity250
## Rugosity750 -0.05858
## Null
## STD.Dpt.1000:STD.Slp.1000 family df logLik AICc delta
## Rugosity500 poisson(log) 5 -2380.883 4771.8 0.00
## Rugosity1000 -0.08218 poisson(log) 5 -2417.702 4845.5 73.64
## Rugosity250 poisson(log) 5 -2430.442 4870.9 99.12
## Rugosity750 poisson(log) 5 -2435.953 4882.0 110.14
## Null poisson(log) 2 -2485.816 4975.6 203.82
## weight
## Rugosity500 1
## Rugosity1000 0
## Rugosity250 0
## Rugosity750 0
## Null 0
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(Rugosity500)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: FCAbundance ~ (1 | Year) + STD.Depth.500 * STD.Slope.500
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 4771.8 4796.7 -2380.9 4761.8 1075
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5311 -0.9536 -0.2128 0.6746 5.6832
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.2325 0.4822
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.15750 0.15421 7.506 6.10e-14 ***
## STD.Depth.500 0.12392 0.02882 4.299 1.71e-05 ***
## STD.Slope.500 -0.34553 0.04262 -8.106 5.22e-16 ***
## STD.Depth.500:STD.Slope.500 -0.05202 0.01983 -2.623 0.0087 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) STD.Dp.500 STD.S.
## STD.Dpt.500 -0.027
## STD.Slp.500 0.070 -0.728
## STD.D.500:S -0.086 0.152 -0.511dist.to.dep.aic <- list("Dist.to.0" = Dist.to.0,
"Dist.to.10" = Dist.to.10,
"Dist.to.20" = Dist.to.20,
"Dist.to.30" = Dist.to.30)
model.sel(dist.to.dep.aic)## Model selection table
## (Int) Dst.to.0m Dst.to.10m Dst.to.20m Dst.to.30m family df
## Dist.to.30 1.123 0.2686 poisson(log) 3
## Dist.to.0 1.125 -0.2583 poisson(log) 3
## Dist.to.20 1.126 0.2517 poisson(log) 3
## Dist.to.10 1.158 0.03017 poisson(log) 3
## logLik AICc delta weight
## Dist.to.30 -2357.706 4721.4 0.00 0.917
## Dist.to.0 -2360.435 4726.9 5.46 0.060
## Dist.to.20 -2361.376 4728.8 7.34 0.023
## Dist.to.10 -2484.016 4974.1 252.62 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'summary(Dist.to.20)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: FCAbundance ~ (1 | Year) + Distance.to.20m
## Data: TurtleData
##
## AIC BIC logLik deviance df.resid
## 4728.8 4743.7 -2361.4 4722.8 1077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3598 -0.9920 -0.1704 0.7198 5.5644
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.2325 0.4822
## Number of obs: 1080, groups: Year, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.12581 0.15361 7.329 2.32e-13 ***
## Distance.to.20m 0.25169 0.01594 15.790 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Dstnc.t.20m -0.027Turt.FC.data1 <- data.frame(FCs = TurtleData$FCAbundance,
Year = as.numeric(TurtleData$Year),
muDep = abs(TurtleData$MEAN.Depth.1000),
muSlop = TurtleData$MEAN.Slope.250,
Rugosity = TurtleData$STD.Slope.500,
dist30 = TurtleData$Distance.to.30m,
shelfW = (TurtleData$Distance.to.20m - TurtleData$Distance.to.10m)
)
ggpairs(Turt.FC.data1, columns = 3:7) + theme_bw()
GlobalModel <- glmer(FCs ~ muDep + muSlop + Rugosity +
dist30 + shelfW + (1|Year), data = Turt.FC.data1, family = poisson)
performance::check_collinearity(GlobalModel)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## muDep 1.08 1.04 0.93
## muSlop 1.74 1.32 0.58
## Rugosity 1.60 1.26 0.63
## dist30 2.66 1.63 0.38
## shelfW 1.52 1.23 0.66GlobalModel <- glmer(FCs ~ muDep + muSlop + Rugosity +
shelfW + (1|Year), data = Turt.FC.data1, family = poisson)
performance::check_collinearity(GlobalModel)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## muDep 1.07 1.04 0.93
## muSlop 1.09 1.04 0.92
## Rugosity 1.31 1.14 0.76
## shelfW 1.25 1.12 0.80dredge(GlobalModel)## Fixed term is "(Intercept)"## Global model call: glmer(formula = FCs ~ muDep + muSlop + Rugosity + shelfW + (1 |
## Year), data = Turt.FC.data1, family = poisson)
## ---
## Model selection table
## (Intrc) muDep muSlp Rgsty shlfW df logLik AICc delta weight
## 16 1.168 -0.06482 -0.1502 -0.1566 0.1285 6 -2323.531 4659.1 0.00 0.657
## 15 1.113 -0.1542 -0.1475 0.1352 5 -2325.194 4660.4 1.30 0.343
## 11 1.121 -0.1804 0.1817 4 -2344.831 4697.7 38.56 0.000
## 12 1.147 -0.02961 -0.1793 0.1796 5 -2344.470 4699.0 39.85 0.000
## 8 1.237 -0.13510 -0.1311 -0.2352 5 -2353.371 4716.8 57.66 0.000
## 7 1.122 -0.1373 -0.2233 4 -2361.025 4730.1 70.95 0.000
## 14 1.217 -0.11230 -0.2147 0.1155 5 -2361.194 4732.4 73.30 0.000
## 13 1.122 -0.1998 0.1272 4 -2366.224 4740.5 81.34 0.000
## 6 1.272 -0.16910 -0.2843 4 -2383.573 4775.2 116.04 0.000
## 5 1.129 -0.2707 3 -2395.585 4797.2 138.05 0.000
## 10 1.197 -0.07077 0.1929 4 -2399.262 4806.6 147.42 0.000
## 9 1.135 0.1979 3 -2401.327 4808.7 149.54 0.000
## 4 1.245 -0.12060 -0.1809 4 -2415.528 4839.1 179.95 0.000
## 3 1.141 -0.1845 3 -2421.862 4849.7 190.60 0.000
## 2 1.296 -0.15990 3 -2474.804 4955.6 296.49 0.000
## 1 1.158 2 -2485.816 4975.6 316.50 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Year'res <- dredge(GlobalModel)## Fixed term is "(Intercept)"importance(res)## muSlop shelfW Rugosity muDep
## Sum of weights: 1.00 1.00 1.00 0.66
## N containing models: 8 8 8 8car::Anova(GlobalModel)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: FCs
## Chisq Df Pr(>Chisq)
## muDep 3.327 1 0.06815 .
## muSlop 73.968 1 < 2.2e-16 ***
## Rugosity 38.080 1 6.791e-10 ***
## shelfW 58.317 1 2.231e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(Turt.FC.data1, aes(muDep, FCs)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.FC.data1, aes(muSlop, FCs)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.FC.data1, aes(Rugosity, FCs)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'
ggplot(Turt.FC.data1, aes(shelfW, FCs)) +
geom_point(size = 1) +
geom_smooth(method = "glm", se = FALSE,
method.args = list(family = "poisson"))## `geom_smooth()` using formula 'y ~ x'